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Abstract 



The problem of close-packed dimers on the honeycomb lattice was solved by 
Kasteleyn in 1963. Here we extend the solution to include interactions between 
neighboring dimers in two spatial lattice directions. The solution is obtained by 
using the method of Bethe ansatz and by converting the dimer problem into a 
five-vertex problem. The complete phase diagram is obtained and it is found 
that a new frozen phase, in which the attracting dimers prevail, arises when the 
interaction is attractive. For repulsive dimer interactions a new first-order line 
separating two frozen phases occurs. The transitions are continuous and the 
critical behavior in the disorder regime is found to be the same as in the case of 
noninteracting dimers characterized by a specific heat exponent a = 1/2. 
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1 Introduction 



An important milestone of the modern theory of lattice statistics is the exact solu- 
tion of the dimer problem obtained by Kasteleyn [jT] and by Fisher . Kasteleyn 
and Fisher considered the problem of close-packed dimers on the simple quartic 
lattice and succeeded in evaluating its generating function in a closed-form ex- 
pression. While the solution shows that close-packed dimers on the square lattice 
do not exhibit a phase transition, Kasteleyn H] later pointed out that dimers on 
the honeycomb lattice do possess phase changes, and that the transitions are 
accompanied by frozen ordered states. The solution, which has since been ana- 
lyzed by one of us [f|, || , can be used to describe domain walls in two dimensions 

In this paper we consider once again close-packed dimers on the honeycomb 
lattice, but now with the introduction of interactions between neighboring dimers 
along two lattice directions. We show that, with the onset of dimer-dimer in- 
teractions, a new ordered phase emerges if the interaction is attractive. For 
repulsive interactions the phase diagram is drastically changed and a tricritical 
point emerges. We deduce locations of all phase boundaries and study its critical 
behavior. 

We analyze interacting dimers by first converting the problem into a five- 
vertex model. For noninteracting dimers this leads to a free-fermion model || 
which can be solved by using the method of Pfaffians Jl|, ]4|, But when interac- 
tions are present the five- vertex model has general vertex weights and the method 
of Pfaffians is no longer applicable. While its solution is in principle obtainable 
from that of the six-vertex model by Sutherland, Yang and Yang announced in 
[§], but details of |9j has not yet been published. Likewise, recently published 
analyses of the general six-vertex model in the regime A < 1 by Nolden |10| and 



in the regime A > 1 by Bukman and Shore |TT| , 12], where A is a parameter 
occurring in the six-vertex model, do not readily translate into the five-vertex 
problem since the five- vertex model corresponds to taking the |A| — > oo limit. 
In fact, it is precisely because of this special situation that the analysis of the 
five-vertex model as a limit of the six-vertex model requires special care. To be 



sure, several authors |7|, [L3| have recently studied the five-vertex model. But 



the five-vertex model considered in [H| is confined to a special regime of the 
parameter space which does not yield the complete complexity of the system. 
The treatment in [[?[], which was aimed to studying domain walls, is more com- 
plete but analyzes the Bethe ansatz solution along a line somewhat different from 
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what we shall present, and is not transparent in extracting relevant information 
on the dimer system. It is therefore useful to have an alternate and self-contained 
analysis of the five-vertex model in the language of dimer statistics. 

We take up this subject matter in the present paper. Our approach is essen- 
tially that of ||, by considering the solution of the Bethe ansatz equations in the 
complex plane. However, we follow the Bethe ansatz solutions closely and explic- 
itly carry out all relevant contour integrations as dictated by relevant physical 
considerations. This leads to a complete and clear picture of the phase diagram 
and critical behavior of the interacting dimer system. Particularly, we find the 
emergence of a new ordered phase for attractive dimer-dimer interactions, and 
the existence of a first-order line terminating at a new kind of tricritical point, 
when the interactions are repulsive. 

The organization of our paper is as follows. The problem of interacting dimers 
is defined in section 1 and mapped into a five- vertex model in section 2. The 
Bethe ansatz equation is set up in section 3, and solved in section 4 in the 
case of noninteracting dimers. In section 5 we analyze the general Bethe ansatz 
equation, obtaining expressions for the free energy and its derivative. This leads 
to the determination of the contour of integrations in section 6 and the complete 
phase diagram in section 7. Finally, the critical behavior is determined in section 
8 by applying perturbation calculations to the free energy. 

2 The five-vertex model 

Consider close-packed dimers on an honeycomb lattice £ which we draw as a 
"brick- wall" shown in Fig. 1. To each dimer along the three edges incident at 
a vertex, one associates a fugacity, or weight, u,v, or w. A vertical u dimer 
and a horizontal v dimer are said to be neighbors if they happen to occupy two 
neighboring sites in the same row. Let two neighboring u and v dimers interact 
with an energy — e and thus possessing a Boltzmann factor 



with A > 1 (A < 1) denoting attractive (repulsive) interactions. Other pairs of 
dimers such as u-u, u-w etc. are not interacting in our model. Then, by replacing 
the two sites inside each dotted box containing a u>-edge in Fig. 1 by a single 
vertex, and regarding a dimer incident to this vertex as a bond covering the 
corresponding lattice edge, the honeycomb lattice C reduces to a simple quartic 
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lattice, and dimer coverings on C lead to vertex configurations of a five-vertex 
model Q. Configurations of the five- vertex model are shown in Fig. 2 in the 
context of a six-vertex model. It is straightforward to verify that we have the 
correspondence 

{u>i, u 2 , u>3, LO4, uj 5 , u e } = {0, w, v, u, V Xuv, V Xuv}. (2.2) 

Here, for definiteness, we assume {u, v, w} > 0. Note that ( |2.2| ) is the most 
general five- vertex model, since one can always take yXuv = y/u 5 ujQ, if u; 5 7^ cu 6 . 
The partition function of a vertex model is defined as 

Z= J2 II u Hv), (2-3) 

config vertices 

where the summation is taken over all allowed vertex configurations, the product 
is over all vertices of the square lattice and lom^ stands for the weight of a vertex 
v . The case of e = or A = 1 leads to the free-fermion model satisfying the 
free-fermion condition ujiU 2 + = uj 5 Uq p. 

For a simple quartic lattice of size M x N, one defines the per-site free energy 

/(w; A)^lim^lnZ, (2.4) 

for the five- vertex model. It follows that the per-site generating function for the 
dimer problem defined in a similar way is //2. 

Ordered States: It is instructive to examine the possible ordered states of the 
dimer lattice. When u, v, or w dominates, the ordered states are those shown 
respectively in Figs. 3a, 3b, and 3c, where the lattice C is completely covered 
by u, v, or w dimers. These are the ordered states occurring in the free-fermion 
case But when A dominates (large positive e), a new ordered state can 
materialize as shown in Fig. 3d. It is this ordering that adds to new features to 
the interacting dimer system. 

3 The Bethe ansatz equation 

To begin with consider the general six-vertex model on a simple quartic lattice 
of M rows and N columns with periodic boundary conditions in both directions. 
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Applying transfer matrix in the vertical direction and using the fact that n, the 
number of empty edges (those not covered by bonds) in a row of vertical edges, is 
conserved, one can evaluate the partition function (|2.3|) using the Bethe ansatz 



141 . The Bethe ansatz formulation for the general six- vertex model has been 



given in ||, |Tf| in a ferroelectric language from which the five-vertex limit does 
not follow straightforwardly. Here, for completeness, we state the Bethe ansatz 



equation for the six-vertex model in terms of the vertex weights JTG 
In the limit of large M, N, one finds 



max [A R (n) + A L (n)] M , (3.1) 

{n} 



with 



= u/i ni — - — — — -) 



AlH = w 4 11 ri — , 3.2 



TV = ( _ 1} n+l 



for n < N, where the n complex numbers Zj, j = 1,2, ...,n are the solution of 
the Bethe ansatz equation 

n(ffefy). •< ' - <«) 

with 

z) = U2UJ4 + ZZ' — (U1U2 + ^3^4 — U 5 Uq)z'. (3.4) 

Note that for fixed 1 < n < N, one has generally r^j different A#(n) and A^(n). 
It is understood that it is the largest ones for each n that are used in ( |3.1| ). We 
remark that a useful parameter occurring in the analysis of the six- vertex model 
is 

_ + ^3^4 ~ ^5^6 

2 v /tUiU;2W3^4 

It is then clear that |A| — ► 00 in the limit of u)\ — > 0, making the five- vertex 
model a very special limit. 

Specializing ( p.2[ ) - (|3~3| ) to the five- vertex model weights fl2.2|) , one obtains 

Afl(n) = (/3w)X;v 

n 

AlH = w^n^i + ^i), (3.5) 
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where the z/s are to be determined from the Bethe ansatz equation 



z 



i=l 



fl(^4r)> J = l,2,...,n (3.6) 



with 

zi = -, x 2 = -A, = -(1 - A). (3.7) 

u u w 

It is clear that A# does not contribute unless n = N. But for n = N the partition 
function can be trivially evaluated. In this case there are no u dimers and hence 
each row of £ is covered completed by v or w dimers and one has 

Z = (w + v ) , 
f(u,v,w;X) = max{ln w, lnt> }. (3.8) 

Alternatively, one can show from (p|) and fl3j) that A R {N) + A L (N) = w N + v N 
from which (|3.8|) also follows. Hence from here on we consider only. 
Combining (|2.4|) and ( |3 . X|) , one has 



f(u,v,w; A) = max/(n), (3.9) 
M 

where 

/(n)= lim ilnArCn). (3.10) 

N^oo iv 

The prescription of the Bethe ansatz is that, for each fixed n, one solves ( |3.6| ) 
for Zj. This leads to generally many sets of solutions. One next picks the set of 
solution which gives the largest f(n) for each n. Then, the free energy ( |3.9|) is 
given by the largest f(n) among all n. We shall refer to the set of Zj which gives 
rise to the final expression of the partition function (|3.9|) the maximal set, and 
the prescription of maximization the maximal principle. 

We first point out some immediate consequences of ( |3.6| ) . First , it is clear that 
if Zj is solution of (|3.6|), then its complex conjugate z* is also a solution^ so that 
the ZjS are distributed symmetrically with respect to the real axis. Secondly, 
for = even at least, the negations of (3 and Zj leave ( |3.6|) unchanged. Thus, 
if Zj is a solution of ( |3.6|) , then — Zj is the solution when (3 is replaced by —(3. 



lr rhis has the consequence that one must write z — \z\e 1 ^ with — 7r < cf> < ir, implying that 
all branch cuts must be taking along the negative real axis. This observation plays a major 
role in ensuing considerations. 
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Combining these two observations, we find the solutions for (3 and —(3 related by 
a simple reflection about the imaginary axis. Finally, multiplying the n equations 
in (|3.6|), one obtains the identity 



n 



N 



1. 



(3.11) 



4 Noninteracting dimers 

It is instructive to apply the Bethe ansatz consideration to the free-fermion case 
of [|J. In this free-fermion case we have {xi,^} = {w/u,v/u}, X — 1, (3 = 0, 
and ( |3.6|) becomes 

zf = (-l) n+1 , j = l,2,...,n. (4.1) 

Thus the Zj's are on the unit circle \z\ = 1, and can take on any n of the N roots 
of (|4.1| ). This gives (^j eigenvalues A L (n) as expected. 

For fixed x±,X2, the maximal set of Zj which gives the largest A^(n) for each 
n is obtained by choosing the n largest \x\ + x%Zj\. Now for any n write 

Zj = e^, a = n/N. (4.2) 

In the limit of large M, N, the Zj's are distributed continuously on the unit circle 
with a uniform angular density N/2tt. For fixed a and write z a = e ian , the n 
largest \xi + x^Zj\ are those given by the Zj's on the arc of the circle extending 
from z* a to z a as shown in Fig. 4. One then obtains from ( |3.5| ) and ( |3.10| ) after 
replacing n by a in the argument 

1 [ z a 

f(a) = lnu + - — r / ln(xi + X2z)dz. (4.3) 

Z7TZ Jz* 



The maximal free energy is therefore, after using fl3.9p , 

max/ 

{a} 

/(«o) 



f(u,v,w;l) = max/(a) 
{«} 



= lnu + — / ln(xi + x 2 e l6 )d9, (4.4) 

where a is determined by 

/'(ao) = hi|xi + x 2 2; ao | = 0. (4.5) 
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If l,xi,X2, or equivalently u,v,w, form a triangle, we have < «o < 1 and 
/ = f(a ) analytic in u,v,w. If l,Xx,x 2 , or equivalently u,v,w, do not form a 
triangle, then there are two possibilities. For x\ + x 2 < 1 or w + v < u, we have 
\x\ + x 2 Zj\ < 1 for all Zj, and as a consequence the maximal set is the empty set 
or, equivalently, a = 0. This leads to 

f(u, v , w; 1) = \nu. u>w + v (4.6) 

For \xi — x 2 \ > 1 or \w — v\ > u, we have \xi + x 2 Zj\ > 1 for all Zj, so that we 
take the maximal set a = 1 and z ao = e ,7r . This leads to after carrying out the 
integration in ( [4.4j ) 

f(u,v,w;l) = lnw + maxln{x2, Xi} 

= max {lnv, In w}. u < \w — v\ (4-7) 

Thus, the phase boundary is 

|xi±X2| = l, or \w±v\=u. (4.8) 



These results are in agreement with |4J. The ordered states ( |4.6|) and (^4.7|) are 
the frozen states shown in Figs. 3a - 3c in which the dimer lattice C is completed 
covered by u, v, or w dimers. The transitions are of second order. We note, in 
particular, that the phase boundary ( |4.8|) is determined by setting f'(ao) = 
at oq = 0, or 7r, the two points where the path of integration in ( |4.3|) either 
just emerges or completes a closed contour. The observation of this mechanism 
underlining the onset of phase transitions proves useful in later considerations. 

5 Analysis of the Bethe ansatz equation 



We now return to the Bethe ansatz equation (|3.6| ). Define a constant C(a, (3) by 

i N 



-ir +i n(i-^), (5.i) 



with ip = argC lying in interval (—tt/N, tc/N]. Then the Bethe ansatz equation 
rU) becomes 

TV 



< = (l-/3^r, (5.2) 



S 



and the iVth root of ( |5.2|) gives a trajectory T on which all solutions Zj must 
reside, 

C^Zj = (1 - /3zj) a e^ } & = 2wj/N, j = 1, 2, N, (5.3) 

where C = C(a,(3) = \C(a,/3)\. The trajectory T is a curve in the complex 
z plane which is symmetric with respect to the real axis and is given by the 
equation 

C\z\ = \l- pz\ a . (r) (5.4) 

However, by diagonalizing the transfer matrix explicitly for N < 18, Noh and 



Kim |17j have found that for (3 > 1 the largest eigenvalue assumes the "bounded 



magnon" ansatz of Noh and Kim |7fl in the form of 



z x = z^ aN - 1 

z i = i = 2,3,---,n, (5.5) 

where \zj\ = 1 for all j. This says that in the thermodynamic limit of iV — ► oo 
and 0,1, one root z\ resides at infinity while all other roots converge to 1/ '/3. 
Using this ansatz, one obtains from ( |5.2|) 



C= 0, 0,1, P > 1, (5.6) 

and from ( |3.5| ) and ( |3.10| ) 



f(a) = (1 - a)\nu + alnv, a ^ 0, 1, ^ > 1. (5.7) 

While we shall make use of ( |5.7| ) to determine the phase boundary, however, to 
make our presentation self-contained we shall proceed for the time being without 
using the anzatz (|5.5[) and the result (|5.7f) . It will be seen that one is led to the 



same phase boundary for (3 > 1. 
Define 

rf = C/|/3|, d c = a a {l-af- a . (5.8) 

Then, by examining solutions of (|5.4|) on the real axis, it is straightforward to 
verify that for (3 > 0, T assumes the topology shown in Figs. 5a - 5c, respectively 
for d > d c ,d = d c , and d < d c . The topology of T for (3 < is deduced from those 
in Fig. 5 by applying a reflection about the imaginary axis, and is shown in Fig. 
6. Note that the constant C can be determined once one point on V is known. 
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Particularly, if T intersects the negative real axis at x 
we have 

C = h(a,R) = ll + P R]a . 

\R\ 



-R, where R > 0, then 
(5.9) 



Note that, since contours can be deformed as long as they do not cross poles and 
move along branch cuts, it is the topology of the contours that is important. 

In the limit of large N, the distribution of zeroes on T becomes continuous. 
Let p(z) be the density of zeroes so that N / p(z)dz over any interval of T gives 
the number of z^s in that interval. Then, using ( |5.3|) one finds 



P(z) 



2m' \z 



a 



(5.10) 



Let r be the segment (s) of V symmetric with respect to the real axis and on 
which the maximal set of z/s resides. For fixed a, we find as in Q4.3|) - (|4.5| ), 



f(a) = \nu + ^-f (- 
2m Jr \ z 



z ^)Hx l+X2 z)d, 



(5.11) 
(5.12) 



f(u,v,w; A) = /(a ), 

where ao is the value of a which maximizes f(ct). 

Note that the a-dependence of the free energy ( 5.11 ) for fixed a now enters 
through both r and p(z). However they obey two constraints. First, the fact 
that there are n z^s implies 



L/ iz)dZ = ^nL {l-J^) dZ 



a. 



In addition, taking the absolute value of Q3.ll ), one obtains 



(5.13) 



(5.14) 



To 



which leads to, in the iV — > oo limit, 



f p{z) In zdz = -—((- 
Jr 2m Jr \z 



a 



-i 



\nzdz = 0. 



(5.15) 



For each fixed a, the two constraints (|5.13 ) and ( 5.15|) together with the maximal 
principle of the free energy are sufficient to determine r . Once r is known, the 
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free energy can be evaluated using ( |5.11| ) and ( p,12| ). In carrying out integrations 
along r , one is aided by the fact that the path of integration can be deformed 
as long as it does not cross poles nor run along branch cuts. Care must be taken, 
however, when T intersects branch cuts. When this happens, integration along 
T can be computed by completing the contour into a closed loop and using the 
Cauchy residue theorem. 

In our discussions below we shall also need to evaluate f'(a). Using (|5.11|) 



for f(a), generally the a dependence comes in through both the path r and 
the explicit dependence of p(z) on a. Let Tq consist of an open path running 
continuously from to z$ and, in addition to the open path, possibly another 
closed contour intersecting the real axis. Since Tq can be freely deformed except 
the terminal points and the intersection points with the branch cut, the derivative 
of ( |5.11 ) with respect to a is derived from three kinds of contributions. First, 



the two terminal points Zo and Zq will move with a. We write 

z = R e ie , y = p(z )(dz /da), (5.16) 

and, due to the fact that the expression (|5.10|) for p(z) contains a factor l/27ri, 
we have 

y* = -p{z*){dz*/da). (5.17) 

Then, the contribution to f'{ot) due to the a-dependence of the terminal points 
gives rise to y ln(xi + x 2 zo) + V* hi(£i + x 2 Zq). Secondly, if r consists of another 
closed contour intersecting the branch cut in ( 5.11| ) at one or two points z r < 0, 



then since the integration path can be deformed, the a-dependence is through 
the intercepts z r only (which now moves along the branch cut), the contribution 
due to z r can be treated as in the above. Namely, we regard two points just above 
and below z r as two terminal points. This leads to a contribution of 2niJ2 r Ur, 
where 

y r = p{z r ){dz T I da) = pure imaginary = — y*, (5.18) 

and this contribution is the same for all branch cuts. Thirdly, there is a contri- 
bution due to the second term in p(z) as shown in ( |5.10|) . Combining the three, 
one obtains 

f'(a) = y\n(x 1 + x 2 z ) +y*ln(x 1 + x 2 z*) 
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where the ± sign is determined by the orientation of To at the intercepts. 

Finally, it is clear from (|5.11|) and ( 5.19 ) that, in carrying out the contour 
integrations for f(a) and /'(a), it is important to determine the location of the 
branch point —x\jx 2 relative to the contour r . First, from the readily verified 
identity 

(1-/3)-^ = !--, (5.20) 

X\ w 

one locates the branch point — X\jx 2 by 

X\jx 2 > (1 — /5) _1 , w > v 

< (1-/TT 1 , w<v, (5.21) 

We remark that from the inequality 

xi/x 2 < (3<0. (5.22) 

one has (see below) for Fig. 6c the inequality R > Ri > \j3\~ 1 > xi/x 2 , implying 
that the points x = —R, —R\ are both on the branch cut. We further relate 
J2rVr, y, and y* by taking the derivatives of (|5.13| ) and (|5.15| ) with respective to 
a, obtaining 



y + y = 

y\nz + y*lnz* + 2iriJ2(±y r ) = 

r 

It now follows from (|5.19|) - ( |5.24|) that we have 



1 + 



dz 



2iri Jr z - 
If \az 



2ni Jr z — (3 1 



P- 1 
dz. 



(5.23) 
(5.24) 



/'(a) =2/ln(— +x 2 )+y*ln(^+x 2 ) + — 



In z — ln(xi + x 2 z) 
z-p- 1 



dz, (5.25) 



provided that y r 's in ( |5.19| ) and ( |5.24j ) are the same, namely, r cuts both branch 
cuts at the same points. In section 7 we shall compute f'(a>) using ( p. 25 ) which 
applies to all cases and all a, including the case that T consists of an open path 
as well as a closed contour. The derivative f'(a) near the phase boundaries will 
also be computed in section 8 by analyzing small perturbations of the free energy. 



12 



6 The contour r at a = 0, 1—, 1, and 1/2 

In our discussions we shall need to evaluate integrals at a = 0, 1—, 1, and 1/2. 
It turns out that the contours r for a = 1 and a = 1— should be considered 
with care. In this section we consider the contour T at these special points. 

(a) a = 0. This is the case that T begins to emerge with very few z/s. The 
constraint ( |5.14|) then dictates that z ~ 1 and, using (|5.4j ) and (|5.14|) , one finds 
C(0,/3) = 1 and T the unit circle. 

(b) a — 1. This is the case of n = N when one picks all iV z/s and hence To = T. 
One can use ( |5.2| ) to obtain 



7iV 



C JV (1, (3)-(-0) 



A' 



z? + ...-l = 0. 



Using (|6.1[) , the constraint (|5.14j) leads to the relation 



N 



or, in the thermodynamic limit, 

C(l,/3) 



One also finds from 



8f) that d c = 1 for a 

d > d c , 

= d r ., 



\0\<1 
\j3\ >1 

1, and hence 

|/3|<1 
|/3|>1. 



(6.1) 



(6.2) 



The trajectory T is, from ( |5.4j ) with C = C(l,/3), 

1 



/3z 



/3 



- 1 



1. 



\(3\<l 



> 1. 



(6.3) 



(6.4) 
(6.5) 



(6.6) 



(6.7) 



The contour ( |6.6| ) for \(3\ < 1 is a circle of radius (1 — (3 2 ) 1 centered at x = 
-f3/(l-[3 2 ) on the real axis, where x = Re(^). Particularly, the circle intersects 
the real axis at x = — R2, R\, where 



R2 



R 1 = (l + P)- 1 , 



(6.8) 
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The contour Q5.7| ) for \(3\ > 1 is the vertical straight line x = 

It is readily verified that by integrating along the contours ( p.6| ) and 
one has 

— ( p(z)\nzdz = 0, \/3\<l (6.9) 

= \p\>l. (6.10) 



Thus, the constraint Q5.13| ) is not satisfied for \(3\ > 1, indicating that the a = 1 



solution is spurious.Q This leads us to consider instead the a = 1— solution as a 
limit to a = 1. 

(c) a = 1 — . 

We consider the cases (3 > and (3 < separately. 

For (3 > 0, the contours are those shown in Fig. 5. In all cases, since the 
branch cuts are on the negative axis, the contour can be deformed to form a 
single closed contour intersecting the real axis at two points and enclosing the 
origin. Let the intersecting point on the negative real axis be x = —R. One finds 
from ( |5.15 ) the relation 

lnh(a,R) = 0, (6.11) 

where a = 1—, h(a,R) is defined in ( ^9|) and is equal to C. Thus, one obtains 
C = 1. In addition, one can solve for R from ( |6.11| ) and obtains 



R = (I-/?)" 1 , (3<l 

R = oo, /3>1. (6.12) 

Thus, the intercept with the negative axis R for /3 < 1 is the same as that given in 
( |6.8| ). Once C = 1 is known, one can compute the location of other intercept (s) 
with the real axis which will be all positive including the intercept (1 + (3)~ l 
given by ( |6.8|) . Note that the second line of (|6.12| ) (and what follows) is what 
one would obtain without making use of the bounded magnon ansatz ( |5.5|) . In 
this case the precise location of the other two intercepts (for f3 > 1, see below) 
does not concern us since the contour can be deformed freely in the x > plane 
as long as it does not cross the pole at (3~ x . 

For (3 < the contours are those shown in Fig. 6. Consider first Fig. 6a 
where the contour intersects the negative real axis at one point at x = —R, one 



2 However, if one carries out the product Ylr z i YVl * ne of (5.14) for finite N, one can 



show that ( 5.14 ) is satisfied 



14 



finds again C = 1 and the two intercepts ( |6.8| ) with R = (1 — In the 

case of Fig. 6c where T intersects the real axis at four points as shown with 
—R < —R\ < —R2 < R3, generally the contour cannot be deformed into a single 
loop due to the presence of the branch cuts. But the integration (|5.15 ) can be 
carried out as in the above yielding 

- In h(a, R x ) + In h(a, R 2 ) + In h(a, R) = 0, (6.13) 

where 

h{a, R) = h{a, R x ) = h{a, R 2 ) = h{a, -R 3 ) = C. (6.14) 
The identities fl6l3D and (jBUD can hold only for C — 1, 

R, = (\(3\ - R 2 = (|/3| + I)" 1 , (6.15) 



which is the same as ( |6.8| ), and both R and -R3 diverging as in the second line of 
( |6.12|) . Note in particular that R 2 coincides with R = (1 — /3) _1 in ( |6.8j ) of Fig. 



5a. Now we have again from (|5.8|) d c = 1 for a = 1—. Thus, in contradistinction 
to (|B~3| ), one finds C(l—,/3) = 1 for all (3 and therefore 

d=j- > d c , \/3\<l (6.16) 

< d c , \p\ > 1. (6.17) 

Thus, the contour T consists of two loops when > 1, with the outer loop 
residing in the infinite regime. 

We note that while ( |6.16| ) is the same as (|6.4|), ( |6.17| ) is different from (|6.5|). 



Finally, one verifies that, by deforming T into a closed contour enclosing the 
origin but not /3 _1 , the constraint ( |5.13| ) is identically satisfied for all {3. 



(d) a = 1/2. 

We are primarily interested in the free energy and its derivatives near a phase 
transition point which, as in the nonintersecting case, occurs when the maximal 
contour either closes or just begin to emerge. Thus, for a = 1/2. we consider 
the cases shown in Figs. 5c and 6c when the contour T consists of two loops and 
T is one of the two loops. But this cannot happen for j3 > 0. When (3 > the 
inner loop and portion of the outer loop are in the x > half plane. Then r in 
the maximal solution of the free energy ( |5.11| ) cannot be either of the two loops, 



since some points on the other loop will have larger values of ln(xi +x 2 z). Hence 
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a = 1/2 can occur only for (3 < 0. This conclusion is also expected on physical 
grounds, that the a = 1/2 ordered state in Fig. 3d dominates only when the 
interaction between neighboring u and v dimers is sufficiently attractive, namely, 
when A is sufficiently large or (3 sufficiently negative. 

For (3 < the contour T intersects the negative real axis (and the branch cut 
of In z) at three points. One verifies after some algebra the identity 



/ p{z)dz = aln|2;o| ^ a, 

J inner Iood 



' inner loop 

for any a. Therefore T must be the outer loop. 

Taking the outer loop as r , we have firstly from the constraint (|5.13| ) 

a = / ( ——\dz = 1 — a, 

2TxiJv a \z z-(3- l J 

leading to the correct value a = 1/2. The constraint ( 5.15Q now yields 

o = -L/ (1_^L_W, 

27rzir V^ z - / 
= \nR-^\n(R + p- 1 ), (6.18) 



or, equivalently, 

This yields the solutions 



R 2 - R - p- 1 = 0. (6.19) 



^'W 1 -^)- (6 - 20) 



Thus, we see that, as expected, (|6.20|) has solution only for sufficiently negative 
(3 < -4. 

Furthermore, when To is a closed contour, one has from (|5.9|) , ( |5.8|) and using 

C(1,P)=\/W\, dc = \, (6-21) 



2" J v " " 2' 
and 

< d c , (3< -4. (6.22) 
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Thus, T consists of two loops only when (3 < —4. In this case V intercepts the 
real axis at the four points — R < —R\ < — R2 < R3 as shown in Fig. 6c and 
determined from C\x\ = J\l — /3x\. This leads to 



R — R+i R-i — R— 




In summary, we have found that, for a = 1— , we have C = 1 and T consists 
of one loop for |/3| < 1 which intersects the real axis at ( |6.8|) , and two loops for 
|/3 1 > 1 intersecting the real axis at 4 points. In the latter case when (5 > 1, 
the contour can be deformed into one loop enclosing the origin and therefore the 
situation is the same as for \/3\ < 1. In the case of (3 < — 1 the outer loop resides 
in the infinite regime. For (5 < 1, the intercept of V on the negative axis and 
closest to the origin is at x = —(1 — /3) -1 . For a = 1/2, one finds C = and 
that, for f3 < —4, T consists of two loops intersecting the real axis at the four 
points given by ( |6.23| ). 

7 The phase diagram 

Ideally one would like to proceed at this point to compute the free energy ( |5.12| ) 
from which the complete thermodynamics of the dimer system can be deter- 
mined. However, as this evaluation involves path integrations which generally 
cannot be put into closed forms, we shall in the next section apply small pertur- 
bations to the free energy near the phase boundaries. We proceed here to first 
determine the phase boundaries and the phase diagram. 

Guided by the analysis of the A = 1 solution of section 4, we expect singular- 
ities of the free energy ( p,12| ) to occur when either the maximal path T contains 
a small emerging segment or when it completes a closed contour. As in the case 
of A = 1, this will happen at a = (segment emerging) and a = 1 (closed 
contour) leading to the ordered states shown in Figs. 3a - 3c. In addition, we 
expect another singularity of the free energy to occur at a® = 1/2 for (3 < —4 
as T completes the outer loop of two closed contours, leading to the order state 
shown in Fig. 3d. In all these cases the free energy f(a) should be the maximal 
solution at a = cto = 0, 1, or 1/2, and the phase boundary is given by /'(ao) = 
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provided that f(ct) is concave at atQ. Once the contours To is known, both the 
free energy ( |5.11| ) and its derivative ( 5.25 ) can be evaluated at a = 0, 1, 1/2. 

Consider first a = 0. This is the case that To begins to emerge at zq = 
RqC 19 ~ 1 and y + y* — 1. It then follows from ( 5.11 ) and ( 5.25|) that we have, 
for /? < 1 at least, 



/(0) 
f(0) 



In w, 



In 



(7.1) 



For a = 1, (or more precisely a = 1 — ) and 1/2, both integrals ( p. 11 ) and 
(|5.25|) can be evaluated using the contours T determined in the preceeding sec- 
tion. We leave details of the evaluations elsewhere, and collect here the results. 

For a = 1 one finds, for (3 > 1, 

/(I) = lBU, 

f'(l) = ln(v/u), 



and, for (3 < 1, 



/(l) = max {In w, In v} 
/'(I) = In' 



In 



■urn — wu(l — A) 
v — w 



Xu 



(7.2) 



w > v 



V > w 



(7.3) 



For a = 1/2 one finds, for /3 < —4, 



1 



2 
In 



ln(Awv) 

vX 
u 



w 



(7.4) 



where R± is given by (|6.20|) . 

To ensure that the free energy f(a>o) is indeed the maximal solution and that 
f'(a ) = is a phase boundary, we need to ascertain that f(a ) is a maximum. 
This turns out to be a delicate matter for a = 0, 1, as detailed calculations show 
that /"(cto) = 0. However, one can proceed as follows. 
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First consider the case (3 > 1. As discussed in Sec. 5 and guided by numerical 
evidence, the ground state is given by the bounded magnon solutions ( |5.5| ) with 
the free energy f(a) given by ( |5.7| ). It follows that the maximal free energy is 

f(u, v , w; A) = max {ln-u, lnt>}, w<v(l — X). (7.5) 

and the phase boundary for (3 > 1 is /(0+) = /(l— ), or 

u = v for w < v(l — A). (7.6) 

Alternately, one can also arrive at the same phase boundary ( |7.6[) without invok- 
ing ( |5.5| ) and ( |5.7| ): Assuming that f{a) is either concave or convex in < a < 1, 
it is easy to see that /'(0) = and /'(l) = cannot be the phase boundary for 
(3 > 1. This follows from the observations that, using ( |7.1| ), 



/'(0) > -> w + Aw > w 

— » v > u (using /3 > 1 or w + Xv < v ) 



and using (FT 



/(I) > 0, (7.7) 



/'(!)< -> w<m 

— > u > w + Xv (again using (3 > 1) 

- /'(0)<0. (7.8) 

Thus, for /3 > 1, in the former case, /'(0) = is not a phase boundary since 
/(0) < /(l), and in the latter case /'(l) = is not a phase boundary since 
/(l) < /(0). Furthermore, relations (|7.7| ) and ( |7.8|) are consistent with the fact 
that the maximal free energy occurs at, respectively, the frozen states ao = 1 and 
cto = 0. It follows that the accompanying transition is between the two frozen 
phases and is thus of first order. 

However, the phase boundary are given correctly by f'(0) = and /'(l) = 
for (3 < 1. This can be seen as follows. 

For /(0) to be a maximum we have always /'(0) < 0, and for /(l) a maximum 
we have /'(l) > 0. We find that, for f3 < 1 or (1 — A)f < w, one has 

/'(0)<0 ^ /'(1)<0 (7.9) 
/'(!)> ^ f(0)>0. (7.10) 
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Thus, under the same convexity and concavity assumption, /(0) and /(l) can 
indeed be the maximum of the free energy. Thus, the phase boundaries are 
/'(0) = and f'{l) = and the accompanying transition is continuous. Explic- 
itly, the phase boundaries are, at a = 0, 

w + Xv = u for (1 — X)v < w (7-11) 

and, at «o = 1, 

A = forw>v (7.12) 

v = w + \u for w + Xv > v > w. (7-13) 



Finally, for the phase boundary at ao = 1/2, we observe from Q7.4Q that, for 
(3 < —4, f'(a) is discontinuous at a = 1/2 with 

f(l-)>f(l+\ ( 7 - 14 ) 



Then f(a) is a maximum at a = 1/2 provided that we have 

/("-)> and f(l+]<0. (7.15) 



The phase boundary is therefore the borderline cases f{\— ) = and f{\+) 
= when f(^) begins to exhibit a maximum. Now, for (3 < —4 and writing 
vX — \(3\w — v, it can be readily verified that we have 

w 

vX- — >v>0. (7.16) 



R 



± 



Thus, the phase boundary at «o = 1/2 is, from ([7.4|) , 



wA- — = VAm^. (7.17) 

Explicitly, ( [7. 1 T| ) can be written as 

uvX 3 - [w 2 + 2w(u + v) + u 2 + v 2 }X 2 + [2w 2 + 2w(u + v ) + uv]X - w 2 = 0, (7.18) 
which reflects its full symmetry with respect to u and v. 
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The phase diagram. Since the vertex weights ( [2.2| ) are arbitrary to an overall 
constant and since there exists an expected u, v symmetry, it is convenient to 
consider the phase diagram in the parameter space u/w,v/w, A. We have found 
the existence of five regimes W (a = 0), U (a — 1), V (a = 1), A (a = 1/2) 
and D (disordered). The regimes U, V, W, A are phases in which dimers are frozen 
in respective configurations of Figs. 3a,b,c,d with values of cto fixed as indicated. 
The phase boundaries are given by ([7.6|) , (|7.11|) , ( [7.1 2|) , (|7.13|) , and (|7.6| ) , leading 
to the phase diagrams shown in Fig. 9. 

For A = 1, the noninteracting case, the phase diagram is shown in Fig. 9a, and 
the boundaries separating regimes W/D, U/D, and V/D are given respectively 

by am, am, am. 

For A > 1 corresponding to attractive interactions between u and v dimers, 
a new ordered phase A with ao = 1/2 arises. A typical phase diagram is shown 
in Fig. 9b, where, in addition to those phase boundaries already present in Fig. 
9a, a new phase boundary ( [7.1 8|) separates regimes D and A. 

For A < 1 corresponding to repulsive interactions between u, v dimers, a 
typical phase diagram is shown in Fig. 9c. In addition to the boundaries already 
present in Fig. 9a, regimes U and V now share a boundary given by ( |7.6| ), namely 
u = v > w/(l — A), across which there is a first-order transition. This leads to 
the existence of a special transition point at u = v = w/(l — A). It is a point 
where two lines of continuous transition merge into a first order line and maybe 
called a kind of tri critical point. However, it is different from ordinary tricritical 
points in that the discontinuity along the first order line does not vanish at that 
point. 



8 The critical behavior and expansions of the 
free energy near phase boundaries 



In this section we derive expansions of the free energy ( 5.11|) for small deviations 
near phase boundaries, and use the expansions to obtain the critical behavior in 
the disorder regime. 

The phase boundaries are characterized by a ~ 0, 1 and 1/2 and, in all cases 
To can be deformed into a a single trajectory extending from a point z$ to its 
complex conjugation zq. Thus, we write 

z = R e ie , z - /T 1 = Ae^ (8.1) 
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where Ro > 0, A > 0, < {9, 0} < ir. We consider 9, ~ or ir, and there are 
three cases to consider. 

(a) a = 0: In this case To is a small arc of radius Ro extending from angle —9 
to 9. 

For < (3Rq < 1, for example, Q5.13 ) can be written as 



a 



9 0-7T 

a . 

7T 7T 



Similarly, for the cases [3Rq > 1 and (3 < 0, we obtain 

a = a — . 

7T 7T 



(8.2) 



(8.3) 



Furthermore, in all cases Rq, A and \f3\ 1 form a triangle implying the relation 
R :A: \f3\~ 1 = sin0 : sin^ : |sin(0-0)|. (8.4) 

Thus, one obtains 

sin 6 , sin 9 



Ro 



(3 sin( 



.4 



(3sm((j)-9) , 



•5) 



where f3 always has the same sign as </> — 9. For given a and /?, either (|8.2| ) or 
( |8.3| ) and the expression of Rq in ( p.5[ ) relates Rq to 9. To determine i?o and 
individually, another relation connecting R Q and 9 is needed. This is provided 
by (|TJ). 

For all cases, (|5.15| ) can be written as 



= In Rn + — 



2tt Jo 



In 



R e^ - (3- 1 



R e ie - (3- 1 



dip. 



(8.6) 



Thus, Ro and 9 can be generally determined, although implicitly. 

Now we specialize the above consideration to small a. When a is small, 9 
is also small. Then, expanding ( |8.1| ) and ( |S.6| ), one obtains R = 1 + 0(9 3 ), 
A= \1- (3' 1 \+0(9 2 ) and 



9 



1 + 



a 



an, 



(8.7) 
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establishing that 9 ~ qltx. 

Using (fT2|), (|3.3|) and the small angle expansion 



, one finds after some 
algebra that in all cases the expansion of the free energy ( |5.11| ) at a = is 



f(a) = \nu + aln(xi + x 2 ) — 



a 3 7r 2 



XiX 2 



+ 0(a 4 



6 V(a;i+a; 2 ) 2 , 

Note that the corresponding expression Eq. (34) of [|7| contains a typographical 
error. 

(b) a = 1: In this case the contour To is almost a closed loop, and can be 
considered as a closed loop T' intercepting the negative real axis at — Ro plus a 
small arc extending from z^ to zq. Now we have always 6 ~ it and, depending 
on whether <fi ~ it {(3 > 0) or (f> ~ {(3 < 0), we have the two cases to consider. 
Thus, ( |5.13|) leads to 



9-tt 



a- 



— 7T 



7T 

6 -TV (j) 

1 H ce- 
il IX 



71 



f3>0 
(3<0 



We again find ( fg.5| ) hold for all j3, and that, ( |5.15| ) leads to 







a 



\nR -\n\f3A\- — 

ITT JO 



In 



Roe* + p- 1 



dip 



(8.9) 



(8.10) 



We now specialize to a = 1— and # = 7r— . Expanding (|8.10|) and (|8.9|) , one 
obtains Ro = (1 - /3)" 1 + 0[(vr - 0) 2 ], A = [(1 — Z?)" 1 + Z?" 1 ) + 0[{it 

a: 



7T - 



1 



+ 1 



a)7r. 



2 ] and 
(8.11) 



It then follows that, after using ( |5.20| ), flO]), ( ft.ll| ) and some lengthy algebra, 
one arrives at the expansion valid for all [3, 

In [x 2 - x x (1 -/?))+ In 

(a - l) 3 7r 2 / x 2 [xi - /3(xi/3 + x 2 )] 



/(«) 



lnt> + (a — 1) 



6 V(l-/5) 2 [x 2 -(l-/3)a;^ 
for _R ~ (1 — > xi/x 2 , 



Inw + (a — 1) In 



;r 2 



(a - 1 ) 3 tt 2 
6~ 



+ 



Xx{x\P + x 2 ) 



[zz-a-^i] 1 



for i? ~ (1 - /?) 1 <a;i/a; 2 . 



(8.12) 



23 



(c) a = 1/2: In this case the contour Tq (after some deformation in the case of 
a = |+) is a closed loop plus a small arc running from z^ to zq, both intercepting 
the negative axis at —Rq for a ~ ~±. Since —R± < (3 < 0, we have always 
0,0 ~ jr. Therefore ( 5.13Q yields 



a + 







TT 



a- 



0. 



TT 



TT 



TT 



One also finds Rq, A given by (|8.5| ). In addition, (|5.15|) now leads to 



= 21ni? -In A- — f 

2tt Jo 



In 



Roe* + (3- 1 



dip. 



Roe**- ) + 13- 1 

Thus, one obtains R = R T + 0[(tt - 9) 2 ], A = R T + /T 1 + 0[(tt - 



1 



a 



2( - -a )tt. 



from which one deduces after some lengthy algebra the expansion 



(8.13) 



(8.14) 



(8.15) 



/(«) 



In (Amu) + ( a 



In 



+ 



/2(a- 1/2) 



3 7T 2 



\3(R T -l/2Y 



where 



D{x u x 2 , (3) = p- 1 + 



x 2 R% 



D(xi,x 2 ,(3), 



x 2 f3 l )x 2 Rl 



(xi - x 2 R T ) 



(xi - x 2 R T ) 2 



[8.16) 



[8.17) 



Here, the upper (lower) sign pertains to a > 1/2 (a < 1/2). This is an extension 
of the corresponding expressions Eqs. (58), (68) and (69) of @. 

The critical behavior. We have obtained expansions of the free energy (|5.11|) 
in the disorder regime near the phase boundaries «o = 0, 1, 1/2 to be given by, 
respectively, (|S.8|) , (|8.12j ), and (|8.16j ). It is now a simple matter to verify that, 
in all cases, the maximal free energy assumes the form 



f[a (t)] = f[a (0)] + c(u,v,w,\)t 3 / 2 , 



(8.18) 



where t is some measure of a small deviation from the phase boundary in the 
parameter space, c(m, v, w, A) is a function regular in t, and a (t) is the value of 
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a determined from the maximal principle. Considered as a vertex model M, for 
example, t can be |T — T c \, where T c is the critical temperature. It then follows 
from (|8.18|) that the transition is of second order (continuous) and the specific 



exponent is a = 1/2. This is exactly the same critical behavior of noninteracting 
dimers ||, ||, and is also the critical behavior expected from the Pokrovsky- 
Talapov type transitions H]. In order to check the internal consistency of our 
results, however, we shall define t > by writing in respective equations for the 
phase boundary w — > w(l +t), or w — > w(l — t) to ensure in the disorder regime. 
We then expect the resulting expression for c(u,v,w, X) to reflect a {u,v} and 
ocq = {0, 1} symmetry. 



To verify ( 8.18|) , we apply the maximal principle to the free energy 



( |8.12| ), and Q . 1 6| ) . Consider first (|8.8|), the expansion of f(a) at a = 0, for 
which the phase boundary is x\ + X2 = 1 or u = w + Xv . Near the phase 
boundary we write w = w(l + 1), where t is small, and determine ao{t) from 

/Wt)] ^ 1 + ^_I^ ( _^).0. (8,9) 



Substituting this ao(t) into (|8.8|) and expanding for small t, one obtains 




f[a (t)}=\n U + ^J^-f\ (8.20) 

o7T V UVX 



This leads to (8.1 



Consider next ( |8.12| ), the expansion of f(a) at a = 1. For the first line of 
j . 1 2|) , the phase boundary is (|7.13|) or v = w + Xu. Near the phase boundary 



we define t by writing w — » w (1 + 1) in ( f7l3|) and obt am 



f'[a (t)} = ln(l-— t) 1 t" ^) ~ X ] 27r V ^2[^i-/3(^i/3 + ^2)] 



This leads to 



Xu ' 2 V(l-/?) 2 [x 2 -(l-/3)x 1 ] 2 , 

0. (8.21} 



/Mt)]=lnu + ^W^f' 2 . (H.22! 
37T V muA 




Note that ( p.20| ) and ( |8.22| ) reflect the expected {u, v} and «o symmetry. 
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For the second line of ( p.l2j ), the phase boundary is x\ — x^jiX — ft) = 1 or 
( |7.12| ). Near the phase boundary we define t by writing w — > w(l — t) in (|7.12|) 
and obtain 



In ( ] 
0. 



2w — u — v 



-t 



w — V 



[a (t) - 1] 2 7T 2 / xi{xxft + x 2 ) 



\[X2 



ft)x x ] 



This leads to 



f[ato(t)] = kiw-t + 



2(2w-u-v) 
37r^ \uv(u + v) 



(8.23) 



^8.24) 



Here the term linear in t comes from the first term in ( |3.24 ) and does not con- 
tribute to the "specific heat" exponent. 

Finally, Consider (|8.16|) , the expansion of the free energy at a = 1/2. Near the 
phase boundary v\—w/R± = V Xuv or ( [7.18] ), we define t by writing w = w(l—t) 
in ( [7.18Q . Then ao(t) is determined from 



f'[<*o{t)] 



Xuv 



2[a (t) - 1/2] 2 t t 2 
(R T - 1/2)3 



D(x h x 2 ,ft) 



= 0, (8.25) 
where D(x\, x%, ft) is given by (|8.17|) . This yields the maximal free energy 

1 



3/4 



f\a Q (t)} = -\n(Xuv) + , . 

J[ VJi 2 v ; 3n\XuvJ \\2(2w + u + v-2Vu^X) 



uv(l - X) 



t 3/2 , (8.26) 



which reflects the proper {u,v} symmetry. Results ( |B.20| ), ( |8.22| ), ( |8.24| ) and 
( |8~26D now confirm ( |57I5D . In writing down flOOP , (ggg) , ( ggg) and (|^|), we 
have used the respective critical conditions to simplify the expressions. 



9 Summary 

We have solved the problem of interacting dimers on the honeycomb lattice 
by solving the equivalent five-vertex model using the method of Bethe ansatz. 
The free energy is given by (|5.11|) and the maximal free energy by Q5.12 ) with 
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r , the contour of integrations, subject to constraints ( |5.13| ) and Q5.15Q . Phase 
transitions are then associated with contours either just emerging or completing a 
closed loop. This leads to the determination of the phase boundaries ( [7.61) , ( |7.11|) , 
(|7.12|) , (|7.13|) and (|7.18|) , and the phase diagrams shown in Fig. 9. We find the 
occurrence of a new frozen ordered phase for attractive dimer interactions, and a 
new first-order line ending at a tricritical point for repulsive dimer interactions. 
We also find, at a — 1, r consist of one loop for \j3\ < 1 and two loops for 
|/3 1 > 1 with the outer loop residing in the infinite regime. But in the latter 
case To can always be deformed into a single loop in computations of the free 
energy and its derivative with respective to a, much simplifying the algebra. At 
a = 1/2, r is found to be the outer loop of two loops, both of which in the 
finite regime, and this occurs only for (3 < —4. We have also evaluated the free 
energy in perturbative expansions near the phase boundary. This leads to the 
determination of the critical behavior in the disorder regime with the specific 
heat exponent a = 1/2. 
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Figure captions 



Fig. 1. The honeycomb lattice drawn as a brick- wall lattice showing relative 
positionings of the u, v, and w dimers. The dotted boxes correspond to vertices 
of the square lattice. 

Fig. 2. The six vertex model and the associated weights. 

Fig. 3. The four possible ordered states, (a) The U phase with n = or a = 
and u dominating, (b) The V phase with n = N or a = 1 and v dominating, (c) 
The W phase with n = N or a = 1 and w dominating, (d) The A phase with 
n = N/2 or a = 1/2 and A dominating. 

Fig. 4. The maximal contour (heavy curve) for noninteracting dimers. 



Fig. 5. Possible contours and the corresponding solutions of (5.4) for (3 > 0. (a) 
d > d c . (b) d = d c . (c) d < d c . 

Fig. 6. Possible contours and the corresponding solutions of ( |5.4| ) for j3 < 0. (a) 

d > d c . (b) d = d c . (c) d < d c . 

Fig. 7. Phase diagrams for fixed A. (a) A = 1, the noninteracting case, (b) 
A > 1, the case of attractive interactions between u, v dimers. A new ordered 
phase A arises for (3 < —4, or v/w > 4/(A — 1). (c) A < 1, the case of repulsive 
interactions between u and v dimers. The U and V regimes share a first-order 
boundary denoted by the heavy line, the circle denotes a tricritical point. 
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